********************************************************************************
*	Do-File for the Replication of the Main Tables, Tables 1-7
*	"Bring a Friend: Leveraging Financial and Peer Support to Improve Women's Reproductive Agency in India"
*	By: S Anukriti, Catalina Herrera-Almanza, Mahesh Karra
*	Date: 17 December 2025
********************************************************************************

***NOTES:
*1. We note that the variable 'solo' throughout our do-files and scripts refers to the Own voucher treatment assignment.

**********************************************
**# B. Data Import and Index Construction
**********************************************

clear all


* Data import
use "$repo/data_endline_final.dta", clear


**## Construction of Z-score Indices for Main Outcomes (Tables 2-5)   

cap drop z_table*

drop if missing(base_school_yrs, base_current_fp, ///
				base_wants_child, base_num_friend, ///
				base_mobi, base_work, base_ever_visit_fp, ///
				base_mod_method, base_phone)
				
				
// Reverse-coding for FP attitudes and pregnancy outcomes to align with treatment effects
cap drop afraid_being_seen_neg
gen afraid_being_seen_neg = 1 - afraid_being_seen
cap drop preg_ever_neg
gen preg_ever_neg = 1 - preg_ever


global outcomes2 ask_peer_adc asked_hmil asked_other asked_sil asked_nonfam
global outcomes3 visit_adc_fp_corr visit_any_fp_new_corr visit_adc_alone /// 
				 visit_adc_w_hhmil_new visit_adc_w_nohmil_new visit_adc_w_sil_new ///
				 visit_adc_w_other_new				 
			 
global outcomes4 heard_mod_method /// 
				 afraid_being_seen_neg  /// 
				 mod_method_el mod_method_since_bl /// 
				 wants_child preg_ever_neg
global outcomes5 outsider_new  /// 
				 atleast1_visithf_VoutHH2 peer_advise_FP_VoHH2
			 
foreach n in 2 3 4 5 {
    
    // Build macro names
    local in_macro outcomes`n'
    local tmp_macro outcomes`n'_tmp
    local zvar z_table`n'

    // Initialize the global for _tmp vars
    global `tmp_macro'

    // Create missflag based on all variables in the group
    cap drop missflag
    egen missflag = rowmiss($`in_macro')

    // Loop over each variable and create _tmp version with conditional missing
    foreach y of global `in_macro' {
        cap drop `y'_tmp
        gen `y'_tmp = `y'
        replace `y'_tmp = . if missflag > 0
        global `tmp_macro' $`tmp_macro' `y'_tmp
    }

    // Create weighted average summary
    cap drop `zvar'
    egen `zvar' = weightave2($`tmp_macro')

    // Clean up
  
}


************************************
**# C. Tables from main paper
************************************

global covbal base_school_yrs base_current_fp base_wants_child base_num_friend
global cov  base_mobi base_work base_ever_visit_fp base_mod_method base_phone

***********************************************************
**## Table 1: Balance at Baseline for the estimation sample 
***********************************************************


* Create global for balance table 1:
global balance base_age base_school_yrs base_w_hindu base_scst base_obc 	///
			   base_ghunghat base_work base_livewMIL base_mobi 				///
			   base_asset_score base_num_friend base_outsider 			    ///
			   base_no_children base_wants_child base_current_fp 			///
			   base_mod_method base_ever_visit_fp base_costFP_MILopposed
 
est drop _all


* Mean . SD of each treatment arm:
eststo control:  estpost  summarize $balance  if treatment == 0
eststo own:  	 estpost  summarize $balance  if treatment == 1 
eststo baf:  	 estpost summarize  $balance  if treatment == 2 

* Pair-differences
eststo diffcontrolown:  estpost ttest $balance if treatment!=2, by(solo) unequal 
* Adding scalars of F-Test of join significance. 
reg solo $balance   if treatment!=2
estadd scalar F_Obs = `e(N)': diffcontrolown
testparm $balance
estadd scalar F_pvalue = r(p): diffcontrolown

eststo diffcontrolbaf:   estpost ttest $balance if treatment!=1, by(baf) unequal 
* Adding scalars of F-Test of join significance. 
reg baf $balance   if treatment!=1
estadd scalar F_Obs = `e(N)': diffcontrolbaf
testparm $balance
estadd scalar F_pvalue = r(p): diffcontrolbaf

eststo diffownbaf:   estpost ttest $balance if treatment!=0, by(baf) unequal 
* Adding scalars of F-Test of join significance. 
reg baf $balance  if treatment!=0
estadd scalar F_Obs = `e(N)': diffownbaf
testparm $balance
estadd scalar F_pvalue = r(p): diffownbaf

* Export Balance Table: 

esttab ///
control own baf  diffcontrolown diffcontrolbaf diffownbaf ///
using "$repo/tables/table1.tex"  , ///
cells("count(pattern(1 1 1 0 0 0 ) fmt(0)) mean(pattern(1 1 1 0 0 0) fmt(3)) b(star pattern(0 0 0 1 1 1) fmt(3))" ". sd(pattern(1 1 1 0 0 0) par fmt(3))" ) ///
label  collabels(none) nonum plain nomtitle noobs not tex starlevels(* 0.1 ** .05  *** .01) ///
stats( F_pvalue F_Obs, label("\midrule F-test of joint significance: p-value" "F-test: Number of observations" ) fmt(%9.3f  %9.0f )) ///
prehead("\begin{tabular}{lccccccccc} \toprule & \multicolumn{2}{c}{Control (C)}   & \multicolumn{2}{c}{Own}  &  \multicolumn{2}{c}{BAF}  & C - Own & C - BAF & Own - BAF \\ \cmidrule(lr){2-3} \cmidrule(lr){4-5} \cmidrule(lr){6-7} \cmidrule(lr){8-8} \cmidrule(lr){9-9} \cmidrule(lr){10-10}  & N & Mean/SD & N & Mean/SD & N & Mean/SD & Diff. & Diff. & Diff. \\ & (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) & (9) \\ \midrule& \multicolumn{2}{c}{\textbf{Control (C)}}   & \multicolumn{2}{c}{\textbf{Own}}  &  \multicolumn{2}{c}{\textbf{BAF}}  & \textbf{C - Own} & \textbf{C - BAF} & \textbf{Own - BAF} \\ \cmidrule(lr){2-3} \cmidrule(lr){4-5} \cmidrule(lr){6-7} \cmidrule(lr){8-8} \cmidrule(lr){9-9} \cmidrule(lr){10-10}  & N & Mean/SD & N & Mean/SD & N & Mean/SD & Diff. & Diff. & Diff. \\ \textbf{Woman's characteristics} & (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) & (9) \\ \midrule") ///
posthead("") ///
postfoot("\bottomrule  \\[-5ex] \end{tabular}") replace



******************************************************************************
**## Table 2: Seeking company for clinic visits during the intervention period 
******************************************************************************

eststo clear
est drop _all

global outcomes ask_peer_adc asked_hmil asked_other asked_sil asked_nonfam ///
				z_table2

foreach y of global outcomes{
eststo r`y'0: reg `y' solo baf $covbal $cov i.base_village_num  , r
test _b[solo]=_b[baf]
	estadd scalar pvalue = r(p)
	sum `y' if treatment == 0 
	estadd scalar controlmean = r(mean)		
}

esttab  rask_peer_adc0   rasked_hmil0  rasked_other0 rasked_sil0 rasked_nonfam0 ///
		rz_table20 ///
using "$repo/tables/table2.tex" , ///
keep(solo baf) ///
label mlabels(none) collabels(none)  se(%9.3f) b(%9.3f)  ///
brackets nonum  noobs nolines  nogaps starlevels( * 0.1 ** .05  *** .01) ///
stats(N  controlmean pvalue , label("\midrule Observations" "Endline control mean" "p-value: Own = BAF" ) fmt(%9.0f %9.3f %9.3f)) ///
prehead("\begin{tabular}{lcccccc} \toprule") ///
posthead(" & \multicolumn{5}{c}{\textbf{Sought company to visit the ADC from:}} &   \\ \cmidrule(lr){2-6} & Someone &  \shortstack{Husband/\\MIL} & \shortstack{Non-husband/\\MIL} & Sister-in-law & Non-relative & \shortstack{Table 2\\index} \\  & (1) & (2) & (3) & (4) & (5) & (6) \\ \midrule") ///
postfoot("\bottomrule \\[-5ex] \end{tabular}") ///
sfmt(%9.2f) varwidth(25)  replace


***************************************************************************************
**## Table 3: clinic Visits for family planning services during the intervention period 
****************************************************************************************


eststo clear
est drop _all

global visits  	visit_adc_fp_corr visit_any_fp_new_corr visit_adc_alone ///
				visit_adc_w_hhmil_new visit_adc_w_nohmil_new visit_adc_w_sil_new ///
				visit_adc_w_other_new ///
				z_table3

foreach y of global visits{
eststo r`y'0: reg `y' solo baf $covbal $cov i.base_village_num  , r

test _b[solo]=_b[baf]
	estadd scalar pvalue = r(p)
	sum `y' if treatment == 0 
	estadd scalar controlmean = r(mean)		
}



esttab  rvisit_adc_fp_corr0 rvisit_any_fp_new_corr0 rvisit_adc_alone0 /// 
		rvisit_adc_w_hhmil_new0 rvisit_adc_w_nohmil_new0 ///
		rvisit_adc_w_sil_new0 rvisit_adc_w_other_new0 ///
		rz_table30 ///
using "$repo/tables/table3.tex" , ///
keep(solo baf) ///
label mlabels(none) collabels(none)  se(%9.3f) b(%9.3f)  ///
brackets nonum  noobs nolines  nogaps starlevels( * 0.1 ** .05  *** .01) ///
stats(N  controlmean pvalue , label("\midrule Observations" "Endline control mean" "p-value: Own = BAF" ) fmt(%9.0f %9.3f %9.3f)) ///
prehead("\begin{tabular}{lcccccccc} \toprule") ///
posthead(" & \multicolumn{2}{c}{} & \multicolumn{5}{c}{\textbf{Visited the ADC:}} &  \\ \cmidrule(lr){4-8} &  \textbf{\shortstack{Visited the\\ ADC}} &  \textbf{\shortstack{Visited any\\clinic}} &   Alone & \shortstack{with husband/\\MIL} & \shortstack{with non-husband/\\MIL} & with SIL & \shortstack{with\\non-relative} & \shortstack{Table 3\\index} \\  & (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) \\ \midrule") ///
postfoot("\bottomrule \\[-5ex] \end{tabular}") ///
sfmt(%9.2f) varwidth(25)  replace


*****************************************************
**## Table 4: Family Planning and fertility outcomes 
*****************************************************

global covbal base_school_yrs base_current_fp base_wants_child base_num_friend
global cov  base_mobi base_work base_ever_visit_fp base_mod_method base_phone


eststo clear
est drop _all

global outcomes heard_mod_method /// // Knowledge
				afraid_being_seen  /// // Attitude
				mod_method_el mod_method_since_bl /// Use of FP
				preg_ever wants_child  /// Fertility
				z_table4

 
foreach y of global outcomes{


eststo r`y'0: reg `y' solo baf $covbal $cov i.base_village_num , r
	test _b[solo]=_b[baf]
	estadd scalar pvalue = r(p)
	sum `y' if treatment == 0 
	estadd scalar controlmean = r(mean)	
}

esttab  _all ///
using "$repo/tables/table4.tex" , ///
keep(solo baf) ///
label mlabels(none) collabels(none)  se(%9.3f) b(%9.3f)  ///
brackets nonum  noobs nolines  nogaps starlevels( * 0.1 ** .05  *** .01) ///
stats(N  controlmean pvalue , label("\midrule Observations" "Endline control mean" "p-value: Own = BAF" ) fmt(%9.0f %9.3f %9.3f)) ///
prehead("\begin{tabular}{lccccccc} \toprule") ///
posthead(" & \multicolumn{1}{c}{\textbf{Knowledge}} &  \multicolumn{1}{c}{\textbf{Attitude}} &  \multicolumn{2}{c}{\textbf{Use of FP}} &  \multicolumn{2}{c}{\textbf{Fertility}} & \\ \cmidrule(lr){2-2} \cmidrule(lr){3-3} \cmidrule(lr){4-5} \cmidrule(lr){6-7}   &  \shortstack{Number of modern\\FP method heard} & \shortstack{Afraid of\\being seen}  & \shortstack{Using a modern\\method at endline}  & \shortstack{Has used a modern\\method since baseline} &  \shortstack{Pregnancy since\\baseline} & \shortstack{Wants another\\child} & \shortstack{Table 4\\index}  \\  & (1) & (2) & (3) & (4) & (5) & (6) & (7)\\ \midrule") ///
postfoot("\bottomrule \\[-5ex] \end{tabular}") ///
sfmt(%9.2f) varwidth(25)  replace


*****************************************************
**## Table 5: Social connections and peer engagement 
*****************************************************

eststo clear
est drop _all	

global outcomes outsider_new  ///
				atleast1_visithf_VoutHH2 peer_advise_FP_VoHH2 ///
				z_table5

foreach y of global outcomes{

eststo r`y': reg `y' solo baf $covbal $cov i.base_village_num  , r
	test _b[solo]=_b[baf]
	estadd scalar pvalue = r(p)
	sum `y' if treatment == 0 
	estadd scalar controlmean = r(mean)	

}

esttab  _all ///
using "$repo/tables/table5.tex" , ///
booktabs fragment nogaps mlabels(none) nonum  nolines collabels(none) label se(%9.3f) b(%9.3f) brackets keep(solo baf) starlevels( * 0.1 ** .05  *** .01) ///
stats(N controlmean pvalue , label("\midrule Observations" "Endline control mean" "p-value: Own = BAF" ) fmt(%9.0f %9.3f  %9.3f)) ///
prehead("\begin{tabular}{lcccc}\toprule") ///
 posthead(" & & \multicolumn{2}{c}{\textbf{Peer engagement}} & \\ \cmidrule(lr){3-4}  & &  \multicolumn{2}{c}{Has close peers outside HH in village that: } \\ \cmidrule(lr){3-4} &  \shortstack{Number of close peers\\outside HH in village} &   \shortstack{Accompanied to\\health facility}  &  \shortstack{Advised woman\\to use FP} &  \shortstack{Table 5\\index}  \\  & (1) & (2) & (3) & (4)  \\  \midrule") ///
 postfoot("\bottomrule \\[-5ex] \end{tabular}") ///
sfmt(%9.2f) varwidth(25)  replace



***********************************************************************************************
**## Table 6: Heterogeneous effects by mother-in-law opposition to family planning at baseline
***********************************************************************************************

eststo clear
est drop _all

global outcomes ask_peer_adc asked_hmil asked_other visit_adc_fp_corr ///
				visit_adc_w_hhmil_new visit_adc_w_nohmil_new ///
				visit_any_fp_new_corr

foreach y of global outcomes{
eststo reg`y': reg `y' 	i.solo##i.base_costFP_MILopposed  ///
						i.baf##i.base_costFP_MILopposed   ///
						c.base_school_yrs#i.base_costFP_MILopposed    ///
						c.base_num_friend#i.base_costFP_MILopposed    ///
						i.base_current_fp##i.base_costFP_MILopposed  ///
						i.base_wants_child##i.base_costFP_MILopposed  ///
						i.base_obc##i.base_costFP_MILopposed  ///
						c.base_mobi#i.base_costFP_MILopposed  ///
						i.base_work##i.base_costFP_MILopposed  ///
						i.base_ever_visit_fp##i.base_costFP_MILopposed   ///
						i.base_phone#i.base_costFP_MILopposed  ///
						i.base_mod_method#i.base_costFP_MILopposed ///
						i.base_village_num#i.base_costFP_MILopposed ,  r
						
test _b[1.solo] = _b[1.baf]
estadd scalar pvalue1 = r(p)

test _b[1.solo#1.base_costFP_MILopposed ] = _b[1.baf#1.base_costFP_MILopposed ]
estadd scalar pvalue2 = r(p)

test _b[1.solo] + _b[1.solo#1.base_costFP_MILopposed ] = 0
estadd scalar pvalue3 = r(p)

test _b[1.baf] + _b[1.baf#1.base_costFP_MILopposed ] = 0
estadd scalar pvalue4 = r(p)

sum `y' if treatment == 0
estadd scalar controlmean = r(mean)	
}


esttab _all ///
using "$repo/tables/table6.tex" , ///
keep(1.solo  1.solo#1.base_costFP_MILopposed  1.baf 1.baf#1.base_costFP_MILopposed   ) ///
order(1.solo  1.baf  1.solo#1.base_costFP_MILopposed 1.baf#1.base_costFP_MILopposed  ) ///
varlabels(1.solo "Own Voucher" 1.baf "BAF Voucher" 1.solo#1.base_costFP_MILopposed "Own Voucher $\times$ MIL opposed to FP" 1.baf#1.base_costFP_MILopposed "BAF Voucher $\times$ MIL opposed to FP") ///
mlabels(none) se(%9.3f) b(%9.3f)  brackets nonum  noobs nolines collabels(none)   nogaps starlevels( * 0.1 ** .05  *** .01)  ///
stats(N  controlmean pvalue1 pvalue2 pvalue3 pvalue4 , label("\midrule Observations" "Endline control mean" "\shortstack[l]{p-values:\\Own = BAF}"   "Own $\times$ MIL opposed = BAF $\times$ MIL opposed" "Own + Own $\times$ MIL opposed = 0"  "BAF + BAF $\times$ MIL opposed = 0"  ) fmt(%9.0f %9.3f %9.3f  %9.3f %9.3f %9.3f)) ///
prehead("\begin{tabular}{lccccccc} \toprule") ///
posthead(" & \multicolumn{3}{c}{\textbf{Sought company to visit the ADC from:}} &  & \multicolumn{2}{c}{\textbf{Visited the ADC with:}} & \\  \cmidrule(lr){2-4} \cmidrule(lr){6-7}   &   Someone & \shortstack{Husband/\\MIL} & \shortstack{Non-husband/\\MIL} & \textbf{\shortstack{Visited\\the ADC}}  & \shortstack{Husband/\\MIL} & \shortstack{Non-husband/\\MIL}  & \textbf{\shortstack{Visited\\any clinic}}  \\  & (1) & (2) & (3) & (4) & (5) & (6) & (7) \\ \midrule") ///
postfoot("\bottomrule \\[-5ex] \end{tabular}") ///
sfmt(%9.2f) varwidth(25)  replace


*****************************************
**## Table 7: Robustness Checks:Inference
*****************************************

global mainoutcomes z_table2 z_table3 z_table4 z_table5


set seed 1234

eststo clear
est drop _all

					
            
foreach y of global mainoutcomes{

* robust
eststo reg_`y': reg `y' solo baf $covbal $cov i.base_village_num  , robust
	test solo
	estadd scalar solo_robust = r(p)
	test baf 
	estadd scalar baf_robust = r(p)
	test _b[solo]=_b[baf]
	estadd scalar pvalue = r(p)

	sum `y' if treatment == 0
	estadd scalar controlmean = r(mean)	

* cluster SE
reg `y' solo baf $covbal $cov i.base_village_num  , cluster(base_village_num)
	test solo 
	estadd scalar solo_cluster = r(p): reg_`y'
	test baf
	estadd scalar baf_cluster = r(p): reg_`y'
	test _b[solo]=_b[baf]
	estadd scalar pvalue_cluster = r(p): reg_`y'

* Wild cluster bootstrapped SE
boottest {solo} {baf}, cluster(base_village_num) nograph  seed(123) 
	estadd scalar solo_wild = r(p_1): reg_`y'
	estadd scalar baf_wild = r(p_2): reg_`y'
	
boottest solo=baf, cluster(base_village_num) nograph  seed(123)
	estadd scalar pvalue_wild = r(p): reg_`y'
}
 
* Computinig MHT using rwolf2 command (Romano-Wolf)

rwolf2 (reg z_table2	solo baf $covbal $cov i.base_village_num, r ) /// 
	   (reg z_table3	solo baf $covbal $cov i.base_village_num, r ) /// 
       (reg z_table4 	solo baf $covbal $cov i.base_village_num, r ) ///  
	   (reg z_table5	solo baf $covbal $cov i.base_village_num, r ), /// 
		indepvars(solo baf , solo baf, solo baf, solo baf) ///
		reps(3000) holm seed(123)

matrix rw_all = e(RW)

* z_table2
scalar solo_mht2 = rw_all[1,3]
scalar baf_mht2 = rw_all[2,3]
estadd scalar solo_mht2 = solo_mht2 : reg_z_table2
estadd scalar baf_mht2 = baf_mht2 	: reg_z_table2

* z_table3 
scalar solo_mht2 = rw_all[3,3]
scalar baf_mht2	 = rw_all[4,3]
estadd scalar solo_mht2 = solo_mht2 : reg_z_table3
estadd scalar baf_mht2 	= baf_mht2 	: reg_z_table3

* z_table4
scalar solo_mht2 = rw_all[5,3]
scalar baf_mht2	 = rw_all[6,3]
estadd scalar solo_mht2 = solo_mht2 : reg_z_table4
estadd scalar baf_mht2  = baf_mht2 	: reg_z_table4

* z_table5
scalar solo_mht2 = rw_all[7,3]
scalar baf_mht2	 = rw_all[8,3]
estadd scalar solo_mht2 = solo_mht2 : reg_z_table5
estadd scalar baf_mht2	= baf_mht2  : reg_z_table5

***** we need to extract the standard errors to calculate the difference between 
***test _b[solo]=_b[baf]
****estadd scalar pvalue = r(p)
* Exporting Table:

esttab 	reg_z_table2 reg_z_table3 reg_z_table4 reg_z_table5 ///
using "$repo/tables/table7_mht.tex" , ///
keep(solo) varlabels(solo "Own Voucher") ///
mlabels(none) b(%9.3f) not nostar nonum  noobs nogaps  eqlabel(none) collabels(none) nolines ///
stats(solo_robust solo_cluster solo_wild solo_mht2, layout(`"(@)"' `"(@)"' `"(@)"' `"(@)"' ) label("\emph{Robust (p-value)}" "\emph{Clustered (p-value)}" "\emph{WC Bootstrap (p-value)}" "\emph{\shortstack{Romano and Wolf (2005a,b)\\MHT Correction (p-value)}}")) ///
prehead("\begin{tabular}{lccccc} \toprule") ///
posthead(" & \multicolumn{4}{c}{\textbf{Z score}} \\ \cmidrule(lr){2-5}  & \textbf{Table 2} & \textbf{Table 3}  & \textbf{Table 4} & \textbf{Table 5} \\  & (1) & (2) & (3) & (4) \\ \midrule")  ///
postfoot("\midrule") ///
 sfmt(%9.2f) varwidth(30) substitute(\_ _)  replace
 
esttab 	reg_z_table2 reg_z_table3 reg_z_table4 reg_z_table5 ///
using "$repo/tables/table7_mht.tex" , ///
keep(baf) varlabels(baf "BAF Voucher") ///
mlabels(none) b(%9.3f) not nostar nonum  noobs nogaps  eqlabel(none) collabels(none) nolines ///
stats(baf_robust baf_cluster baf_wild baf_mht2 N controlmean pvalue pvalue_cluster pvalue_wild, layout(`"(@)"' `"(@)"' `"(@)"' `"(@)"' `"@"' `"@"' `"@"' `"@"' `"@"' ) label("\emph{Robust (p-value)}" "\emph{Clustered (p-value)}" "\emph{WC Bootstrap (p-value)}" "\emph{\shortstack{Romano and Wolf (2005a,b)\\MHT Correction (p-value)}}" "\midrule Observations" "Endline control mean" "\emph{Robust (p-value: Own = BAF)}" "\emph{Clustered (p-value: Own = BAF)}" "\emph{WC Bootstrap (p-value: Own = BAF)}")) ///
prehead("") ///
posthead(" ")  ///
postfoot("\bottomrule \\[-5ex] \end{tabular}") ///
sfmt(%9.2f) varwidth(30)  substitute(\_ _)  append


* Anderson q-values

preserve

// SOLO
* Initialize matrix
local nvars : word count $mainoutcomes
matrix y_s = J(`nvars', 1, .)

* Loop through outcome variables
local i = 1
foreach var of global mainoutcomes {
    quietly reg `var' solo baf $covbal $cov i.base_village_num, r
    test solo
    matrix y_s[`i', 1] = r(p)
    local ++i
}

* Output matrix to CSV
matrix list y_s
mat2txt, matrix(y_s) saving("$repo/tables/table7_pvalue_solo.csv") replace

***********************Post-estimation command using the standard errors of the Anderson-q values 
test _b[solo]=_b[baf]
estadd scalar pvalue = r(p)

* Save matrix as .dta
drop _all
svmat double y_s
rename y_s1 pval
save "$repo/tables/table7_pvalue_solo.dta", replace

restore

// BAF
preserve

* Initialize matrix
local nvars : word count $mainoutcomes
matrix y_b = J(`nvars', 1, .)

* Loop through outcome variables
local i = 1
foreach var of global mainoutcomes {
    quietly reg `var' solo baf $covbal $cov i.base_village_num, r
    test baf
    matrix y_b[`i', 1] = r(p)
    local ++i
}

* Output matrix to CSV
matrix list y_b
mat2txt, matrix(y_b) saving("$repo/tables/table7_pvalue_baf.csv") replace

* Save matrix as .dta
drop _all
svmat double y_b
rename y_b1 pval
save "$repo/tables/table7_pvalue_baf.dta", replace

restore


// Joint BAF = OWN
preserve

* Initialize matrix
local nvars : word count $mainoutcomes
matrix y_b = J(`nvars', 1, .)

* Loop through outcome variables
local i = 1
foreach var of global mainoutcomes {
    quietly reg `var' solo baf $covbal $cov i.base_village_num, r
    test _b[solo]=_b[baf]
    matrix y_b[`i', 1] = r(p)
    local ++i
}

* Output matrix to CSV
matrix list y_b
mat2txt, matrix(y_b) saving("$repo/tables/table7_pvalue_joint.csv") replace

* Save matrix as .dta
drop _all
svmat double y_b
rename y_b1 pval
save "$repo/tables/table7_pvalue_joint.dta", replace

restore

**** SOLO: Now use Michael Anderson's code for sharpened q-values
preserve

use "$repo/tables/table7_pvalue_solo.dta", clear
set more off

* Collect the total number of p-values tested

quietly sum pval
local totalpvals = r(N)

* Sort the p-values in ascending order and generate a variable that codes each p-value's rank

quietly gen int original_sorting_order = _n
quietly sort pval
quietly gen int rank = _n if pval~=.

* Set the initial counter to 1 

local qval = 1

* Generate the variable that will contain the BKY (2006) sharpened q-values

gen bky06_qval = 1 if pval~=.

* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses are rejected at q = 0.998, etc.  The loop ends by checking which hypotheses are rejected at q = 0.001.


while `qval' > 0 {
	* First Stage
	* Generate the adjusted first stage q level we are testing: q' = q/1+q
	local qval_adj = `qval'/(1+`qval')
	* Generate value q'*r/M
	gen fdr_temp1 = `qval_adj'*rank/`totalpvals'
	* Generate binary variable checking condition p(r) <= q'*r/M
	gen reject_temp1 = (fdr_temp1>=pval) if pval~=.
	* Generate variable containing p-value ranks for all p-values that meet above condition
	gen reject_rank1 = reject_temp1*rank
	* Record the rank of the largest p-value that meets above condition
	egen total_rejected1 = max(reject_rank1)

	* Second Stage
	* Generate the second stage q level that accounts for hypotheses rejected in first stage: q_2st = q'*(M/m0)
	local qval_2st = `qval_adj'*(`totalpvals'/(`totalpvals'-total_rejected1[1]))
	* Generate value q_2st*r/M
	gen fdr_temp2 = `qval_2st'*rank/`totalpvals'
	* Generate binary variable checking condition p(r) <= q_2st*r/M
	gen reject_temp2 = (fdr_temp2>=pval) if pval~=.
	* Generate variable containing p-value ranks for all p-values that meet above condition
	gen reject_rank2 = reject_temp2*rank
	* Record the rank of the largest p-value that meets above condition
	egen total_rejected2 = max(reject_rank2)

	* A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition
	replace bky06_qval = `qval' if rank <= total_rejected2 & rank~=.
	* Reduce q by 0.001 and repeat loop
	drop fdr_temp* reject_temp* reject_rank* total_rejected*
	local qval = `qval' - .001
}
	

quietly sort original_sorting_order
pause off

display "Code has completed."
display "Benjamini Krieger Yekutieli (2006) sharpened q-vals are in variable 'bky06_qval'"
display	"Sorting order is the same as the original vector of p-values"

keep pval original_sorting_order rank bky06_qval
save "$repo/tables/table7_sharpenedqvals_solo.dta", replace

restore


**** BAF: Now use Michael Anderson's code for sharpened q-values
preserve

use "$repo/tables/table7_pvalue_baf.dta", clear
set more off

* Collect the total number of p-values tested

quietly sum pval
local totalpvals = r(N)

* Sort the p-values in ascending order and generate a variable that codes each p-value's rank

quietly gen int original_sorting_order = _n
quietly sort pval
quietly gen int rank = _n if pval~=.

* Set the initial counter to 1 

local qval = 1

* Generate the variable that will contain the BKY (2006) sharpened q-values

gen bky06_qval = 1 if pval~=.

* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses are rejected at q = 0.998, etc.  The loop ends by checking which hypotheses are rejected at q = 0.001.


while `qval' > 0 {
	* First Stage
	* Generate the adjusted first stage q level we are testing: q' = q/1+q
	local qval_adj = `qval'/(1+`qval')
	* Generate value q'*r/M
	gen fdr_temp1 = `qval_adj'*rank/`totalpvals'
	* Generate binary variable checking condition p(r) <= q'*r/M
	gen reject_temp1 = (fdr_temp1>=pval) if pval~=.
	* Generate variable containing p-value ranks for all p-values that meet above condition
	gen reject_rank1 = reject_temp1*rank
	* Record the rank of the largest p-value that meets above condition
	egen total_rejected1 = max(reject_rank1)

	* Second Stage
	* Generate the second stage q level that accounts for hypotheses rejected in first stage: q_2st = q'*(M/m0)
	local qval_2st = `qval_adj'*(`totalpvals'/(`totalpvals'-total_rejected1[1]))
	* Generate value q_2st*r/M
	gen fdr_temp2 = `qval_2st'*rank/`totalpvals'
	* Generate binary variable checking condition p(r) <= q_2st*r/M
	gen reject_temp2 = (fdr_temp2>=pval) if pval~=.
	* Generate variable containing p-value ranks for all p-values that meet above condition
	gen reject_rank2 = reject_temp2*rank
	* Record the rank of the largest p-value that meets above condition
	egen total_rejected2 = max(reject_rank2)

	* A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition
	replace bky06_qval = `qval' if rank <= total_rejected2 & rank~=.
	* Reduce q by 0.001 and repeat loop
	drop fdr_temp* reject_temp* reject_rank* total_rejected*
	local qval = `qval' - .001
}
	

quietly sort original_sorting_order
pause off

display "Code has completed."
display "Benjamini Krieger Yekutieli (2006) sharpened q-vals are in variable 'bky06_qval'"
display	"Sorting order is the same as the original vector of p-values"

keep pval original_sorting_order rank bky06_qval
save "$repo/tables/table7_sharpenedqvals_baf.dta", replace

restore


**** BAF = OWN: Now use Michael Anderson's code for sharpened q-values
preserve

use "$repo/tables/table7_pvalue_joint.dta", clear
set more off

* Collect the total number of p-values tested

quietly sum pval
local totalpvals = r(N)

* Sort the p-values in ascending order and generate a variable that codes each p-value's rank

quietly gen int original_sorting_order = _n
quietly sort pval
quietly gen int rank = _n if pval~=.

* Set the initial counter to 1 

local qval = 1

* Generate the variable that will contain the BKY (2006) sharpened q-values

gen bky06_qval = 1 if pval~=.

* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses are rejected at q = 0.998, etc.  The loop ends by checking which hypotheses are rejected at q = 0.001.


while `qval' > 0 {
	* First Stage
	* Generate the adjusted first stage q level we are testing: q' = q/1+q
	local qval_adj = `qval'/(1+`qval')
	* Generate value q'*r/M
	gen fdr_temp1 = `qval_adj'*rank/`totalpvals'
	* Generate binary variable checking condition p(r) <= q'*r/M
	gen reject_temp1 = (fdr_temp1>=pval) if pval~=.
	* Generate variable containing p-value ranks for all p-values that meet above condition
	gen reject_rank1 = reject_temp1*rank
	* Record the rank of the largest p-value that meets above condition
	egen total_rejected1 = max(reject_rank1)

	* Second Stage
	* Generate the second stage q level that accounts for hypotheses rejected in first stage: q_2st = q'*(M/m0)
	local qval_2st = `qval_adj'*(`totalpvals'/(`totalpvals'-total_rejected1[1]))
	* Generate value q_2st*r/M
	gen fdr_temp2 = `qval_2st'*rank/`totalpvals'
	* Generate binary variable checking condition p(r) <= q_2st*r/M
	gen reject_temp2 = (fdr_temp2>=pval) if pval~=.
	* Generate variable containing p-value ranks for all p-values that meet above condition
	gen reject_rank2 = reject_temp2*rank
	* Record the rank of the largest p-value that meets above condition
	egen total_rejected2 = max(reject_rank2)

	* A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition
	replace bky06_qval = `qval' if rank <= total_rejected2 & rank~=.
	* Reduce q by 0.001 and repeat loop
	drop fdr_temp* reject_temp* reject_rank* total_rejected*
	local qval = `qval' - .001
}
	

quietly sort original_sorting_order
pause off

display "Code has completed."
display "Benjamini Krieger Yekutieli (2006) sharpened q-vals are in variable 'bky06_qval'"
display	"Sorting order is the same as the original vector of p-values"

keep pval original_sorting_order rank bky06_qval
save "$repo/tables/table7_sharpenedqvals_joint.dta", replace

restore

* Please refer to the R script file to construct final formatted table7.tex file.

// Clear everything

clear all